Linear transport equations
The simplest case of the one-dimensional version of the transport equation,
is when the velocity field \(a\) is constant. The Cauchy problem is solved using the method of characteristics, where the solution is constant along the characteristic curves \(x = x_0 + a t\). The solution is then
for any \((x,t) \in \R \times \R_+\). We see that the initial data is transported with the velocity \(a\). For more general cases, the characteristics may not be possible to solve explicitly. However, we can obtain some information of the sulutions with the following a priori energy estimate:
Lemma
Assume \(U\) is a smooth solution of the transport equation decaying to zero at infinity for all \(t \in R_+\) and \(a \in C^1(\R, \R_+)\). Then, \(U\) satisfies the energy bound
Proof
Follows by multiplying the transport equation by \(U\) and integrating over space. Then, use that \(U\) decays to zero at infinity and apply Grönwall's inequality.
The lemma shows that the energy is bounded. Using another functional, the assumptions on \(U\) can be relaxed:
Lemma
Assume \(U\) is a smooth bounded solution of the transport equation. Then, we have
for any \(t > 0\).
Proof
For any \((x, t) \in \R \times \R_+\), there exists \(\xi \in \R\) such that \(U(x, t) = U_0(\xi)\)
Finite difference schemes for the transport equation
For some velocity fields \(a(x,t)\), it may not be possible to derive an explicit formula for the characteristic equation. We thus use numerical methods to approximate the solutions of the transport equation.
Discretization
For simplicity, assume that the velocity field is positive. AS \(\R\) is unbounded, we truncate the domain into \(\Omega = [x_L, X_R]\). Thus, we must impose boundary conditions, which will be discussed below. For simplicity, we use a uniform mesh of mesh size \(\Delta x\) with \(N+1\) points \(x_j\):
We further choose some terminal time \(T\) and divide into \(M+1\) points \(t^n = n \Delta t\). We set the initial approximaition \(U_j^0 := U_0(x_j)\) and update the next approximation \(U_j^{n+1}\) using a finite difference scheme.
Centered finite difference scheme
One such scheme is the forward difference in time and cetral difference in space approximating \((\ref{eq:transport_const})\):
Example
We consider the domain \([0, 1]\) with initial data
u0(x) = sin(2*pi*x)
and \(a = 1\). Since the data is periodic, we impose periodic boundary conditions. Numerically, we implement this by setting
Thus, on the boundary, \(j = 0, N\), we have
For the first time step \(n = 1\), we have
Using a grid of \(50\) mesh points, simulating to time \(T = 3\), we get the following result:
x_L, x_R = 0, 1
T = 3
N = 100
dx = (x_R - x_L)/N
dt = 1 * dx
x = x_L+dx:dx:x_R
t = dt:dt:T
U = central_difference.central_difference_scheme(x, t, u0)
Viz.animate_solution(U,
(x, t) -> u0(x-t),
x, t)
After some time, the solution diverges. Intuitively, the information should propagate from left to right. However, the scheme uses information from both sides. This can be explained rigorously using the the discrete energy
We know that the exact solution have bounded energy. We say that a scheme is energy stable if \(E^n \leq E^0\) for all \(n\).
Lemma
Let \(U_j^n\) be the solutions computed with the central difference scheme. Then, we have the following dicsrete energy estimate:
Thus, the energy grows at each time step for any choice of \(\Delta x\) and \(\Delta t\).
Proof
Similar to the proof of the continuous energy estimate and using the identity
Upwind scheme
To respect the flow of information, we can use forward and backward differences in space depending on the direction of propagation of information, i.e.
Information is "carried by the wind", hence the name upwind. The above equations can be written as
We see we have the central difference scheme with a diffusion term; the right hand side approximates \(\frac{\Delta x \abs{a}}{2} U_xx\). Thus, it adds numerical viscosity to the unstable central difference scheme, which will play a crucial role later.
Now, we do the same numerical experiment as previously with \(a = 1\):
x_L, x_R = 0, 1
T = 1
N = 50
dx = (x_R - x_L)/N
dt = 1.3 * dx
x = x_L + dx:dx:x_R
t = dt:dt:T
a(x, t) = 1
U = upwind.upwind_scheme(x, t, a, u0)
Viz.animate_solution(U,
(x, t) -> u0(x-t),
x, t)
x_L, x_R = 0, 1
T = 1
N = 100
dx = (x_R - x_L)/N
dt = 0.7 * dx
x = x_L+dx:dx:x_R
t = dt:dt:T
a(x, t) = 1
U = upwind.upwind_scheme(x, t, a, u0)
Viz.animate_solution(U,
(x, t) -> u0(x-t),
x, t)
We see that stability depends on the relation \(\frac{\Delta t}{\Delta x}\).
Lemma
If the mesh parameters satisfy the condition
then the upwind solution satisfies the estimate
so the scheme is conditionally stable.
Proof
Start similarly to the proof of the unconditional unstability of the central difference scheme and use the mentioned identity several times.
The above condition is called A CFL condition. We also have \(L^1\) and \(L^\infty\) stability:
Lemma
Assume the above CFL condition holds. Then, the solutions of the upwind scheme satisfy
Proof
The first inequality follows from \(U_j^{n+1}\)being a convex combination of \(U_{j-1}^n\), \(U_j^n\) and \(U_{j+1}^n\).
Discontinuous initial data
We now consider the transport equation with \(a = 1\) in the domain \([0, 1]\) with initial data
u0(x) = x .< 0.5 ? 2 : 1
a(x, t) = 1
This yields the discontinuous solution \(U(x, t) = U_0(x - t)\). We use the upwind scheme with \(N = 50\) and \(N = 200\) mesh points.
x_L, x_R = 0, 1
T = 0.25
N = 50
dx = (x_R - x_L)/N
dt = 0.9 * dx
x = x_L+dx:dx:x_R
t = dt:dt:T
U = upwind.upwind_scheme(x, t, a, u0)
Viz.animate_solution(U,
(x, t) -> u0(x-t),
x, t)
x_L, x_R = 0, 1
T = 0.25
N = 200
dx = (x_R - x_L)/N
dt = 0.9 * dx
x = x_L+dx:dx:x_R
t = dt:dt:T
U = upwind.upwind_scheme(x, t, a, u0)
Viz.animate_solution(U,
(x, t) -> u0(x-t),
x, t)